######################################################
#### Replication materials for Schub ISQ
#### "When Prospective Leader Turnover Promotes Peace"
######################################################

## Outline
  # I.   Manuscript Figures 2-3
  # II.  Online Appendix Table A6 and Figure A4
  # III. Online Appendix Figures A3 and A5-A7

###############################
# Set working directory and load data
library(foreign)

setwd("/Users/robertschub/Dropbox/Rob/RESEARCH IDEAS/10. KickingTheCanFutureUncertainty/ISQ_Final/Replication")

datars<-read.dta("MainData_Schub_ISQ_LeaderTurnover_v2.dta")

###############################
# I. Manuscript Figures 2-3

##########
## Manuscript Figure 2 Left Panel 
##########

par(mar=c(3,2,1,1))
par(oma=c(2,1,0,0))
plot(jitter(datars$sol2y10[datars$milgdpest1>10 & datars$war==0]-0.2,factor=0.5), jitter(datars$warinitstartnxtyr[datars$milgdpest1>10 & datars$war==0],factor=0.2),col="grey50",pch=19,xlim=c(0,16),xaxt="n",yaxt="n",ylab="",xlab="")
points(jitter(datars$sol2y10[datars$milgdpest1<=10 & datars$war==0]+0.2,factor=0.5), jitter(datars$warinitstartnxtyr[datars$milgdpest1<=10 & datars$war==0],factor=0.2),col="black",pch=9)
text(6,0.4,"low military spending",cex=1.3)
text(6,0.6,"high military spending",col="grey35",cex=1.3)
points(2.5,0.6,pch=19,cex=2,col="grey50")
points(2.5,0.4,pch=9,cex=2)
axis(1,at=seq(0,20,by=5),cex.axis=1.3)
axis(2,at=c(0,1),labels=c("peace","war"),cex.axis=1.6)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)



##########
## Manuscript Figure 2 Right Panel 
##########

  ## Mean war occurence when spending is HIGH with variation in turnover
h0<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1>10 & datars$war==0 & datars$triturnover==0]))
h1<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1>10 & datars$war==0 & datars$triturnover==1]))
h2<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1>10 & datars$war==0 & datars$triturnover==2]))
  ## Mean war occurence when spending is LOW with variation in turnover
l0<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1<=10 & datars$war==0 & datars$triturnover==0]))
l1<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1<=10 & datars$war==0 & datars$triturnover==1]))
l2<-mean(na.omit(datars$warinitstartnxtyr[datars$milgdpest1<=10 & datars$war==0 & datars$triturnover==2]))

turns<-c(0,1,2)
high<-c(h0,h1,h2)
low<-c(l0,l1,l2)

par(mar=c(3,4,1,2))
par(oma=c(2,1,0,0))
plot(turns,high,pch=19,col="grey50",cex=3.8,ylim=c(0,0.025),yaxt="n",ylab="",xaxt="n",xlab="")
points(turns,low,pch=9,col="black",cex=3.8)
axis(2,at=seq(0,0.025,by=0.005),labels=c(0,0.5,1,1.5,2,2.5),cex.axis=1.4)
mtext("War Occurrence (%)",side=2,outer=F,line=2.6,cex=1.7)
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)
text(0.42,0.022,"high military spending",col="grey35",cex=1.3)
text(0.42,0.006,"low military spending",col="black",cex=1.3)
text(0.12,0.025,"n=372",cex=1.1)
text(0.1,0.0097,"2742",cex=1.1)
text(1.1,0.0184,"154",cex=1.1)
text(1.1,0.0067,"1047",cex=1.1)
text(1.9,0.0064,"1288",cex=1.1)
text(1.9,0.0008,"124",cex=1.1)


##########
## Manuscript Figure 3 
##########

  ## Quantities based on Table 1 Model 6: see "Schub_ISQ_LeaderTurnover_MainAnalysis.do"
fd<-c(-0.204,-1.434)
lb<-c(-0.549,-3.045)
ub<-c(0.248,-0.450)
sec<-c(1:2)

par(oma=c(1,9,1,1))
par(mar=c(3,3,3,3))
plot(fd,sec,xlim=c(-4.2,1.1),ylim=c(0.5,2.5),col="white",xaxt="n",xlab="",ylab="",yaxt="n")
segments(lb,sec,ub,sec,lwd=3,col="grey50")
points(fd,sec,pch=20,cex=2)
abline(v=0,lwd=2,lty=3)
axis(1,at=c(-4:3),cex.axis=1.4)
mtext("Marginal Effect of Previous Leader Turnover",3,outer=T,line=-2,cex=1.6)
axis(2,at=c(1,2),labels=c("Low Mil Spending","High Mil Spending"),cex.axis=1.4,las=2)
mtext("Percentage Point Change",1,outer=F,line=2.4,cex=1.4)




###############################
# II. Online Appendix Table A6 and Figure A4

##########
## Table A6
##########

  ## Penalized (Firth) Logit
library(logistf)
library(brglm)
library(arm)
library(stargazer)

  ## Models 1-6
m1<-brglm(warinitstartnxtyr ~ triturnover,data=datars[datars$war==0 & datars$milgdpest1>10,],family="binomial",method="brglm.fit")
m2<-brglm(warinitstartnxtyr ~ triturnover +  relcap + cont + alliancedummy + pceyrs + pceyrs2 + pceyrs3,data=datars[datars$war==0 & datars$milgdpest1>10,],family="binomial",method="brglm.fit")
m3<-brglm(warinitstartnxtyr ~ triturnover,data=datars[datars$war==0 & datars$milgdpest1<=10,],family="binomial",method="brglm.fit")
m4<-brglm(warinitstartnxtyr ~ triturnover +  relcap + cont + alliancedummy + pceyrs + pceyrs2 + pceyrs3,data=datars[datars$war==0 & datars$milgdpest1<=10,],family="binomial",method="brglm.fit")
m5<- brglm(warinitstartnxtyr ~ triturnover + milgdpest1binary10 + triturnover:milgdpest1binary10 +  relcap + cont + alliancedummy + pceyrs + pceyrs2 + pceyrs3,data=datars[datars$war==0,],family="binomial",method="brglm.fit")
base <- brglm(warinitstartnxtyr ~ triturnover + milgdpest1binary10 + triturnover:milgdpest1binary10 +  relcap + cont + alliancedummy + new_polity22 + pceyrs + pceyrs2 + pceyrs3,data=datars[datars$war==0,],family="binomial",method="brglm.fit")

stargazer(m1,m2,m3,m4,m5,base, align=F,style="apsr",omit.stat=c("f","ser","adj.rsq","rsq"),digits=2)


##########
## Figure A4 based on Model 6 above
##########

beta.draw<-mvrnorm(1000,mu=coef(base),Sigma=vcov(base))
dat<-cbind(1,base$model[,-1])
cov.i<-apply(dat,MARGIN=2,FUN=median)
cov<-c(cov.i,0)

  # HIGH mil spend
cov.on<-cov
cov.on[3]<-1 #set mil spend high
cov.onh<-cov.on
cov.onl<-cov.on
cov.onh[2]<-2 #set turnover to 2
cov.onh[11]<-2 #set interaction to 2
cov.onl[2]<-0 #set turnover to 0
cov.onl[11]<-0 #set interaction to 0
cov.onh
cov.onl
fd.on<-plogis(cov.onh%*%t(beta.draw))-plogis(cov.onl%*%t(beta.draw))
meon<-mean(fd.on)
meonl<-quantile(fd.on,0.025)
meonh<-quantile(fd.on,0.975)

  # LOW mil spend
cov.off<-cov
cov.off[3]<-0 #set mil spend low
cov.offh<-cov.off
cov.offl<-cov.off
cov.offh[2]<-2 #set turnover to 2
cov.offh[11]<-0 #set interaction to 0
cov.offl[2]<-0 #set turnover to 0
cov.offl[11]<-0 #set interaction to 0
cov.offh
cov.offl
fd.off<-plogis(cov.offh%*%t(beta.draw))-plogis(cov.offl%*%t(beta.draw))
meoff<-mean(fd.off)
meoffl<-quantile(fd.off,0.025)
meoffh<-quantile(fd.off,0.975)

me<-c(meon,meoff)
mel<-c(meonl,meoffl)
meh<-c(meonh,meoffh)
ct<-c(2,1)

par(oma=c(1,9,1,1))
par(mar=c(3,3,3,3))
plot(me,ct,xlim=c(-0.04,0.01),ylim=c(0.8,2.2),pch=19,cex=2,yaxt="n",ylab="",xlab="")
segments(mel,ct,meh,ct,lwd=3,col="grey50")
points(me,ct,pch=20,cex=2)
abline(v=0,lwd=2,lty=3)
mtext("Penalized Logit Marginal Effects",side=3,outer=T,line=-2,cex=1.6)
axis(2,at=c(1,2),labels=c("Low Mil Spending","High Mil Spending"),cex.axis=1.4,las=2)
mtext("Percentage Point Change",1,outer=F,line=2.4,cex=1.4)


###############################
# III. Online Appendix Figures A3 and A5-A7

##########
## Figure A3 
##########

  ## Quantities based on Table A3 Model 6: see "Schub_ISQ_LeaderTurnover_MainAnalysis.do"
  
fd<-c(-0.204,-1.128)
lb<-c(-0.567,-2.728)
ub<-c(0.212,-0.158)
sec<-c(1:2)

par(oma=c(1,9,1,1))
par(mar=c(3,3,3,3))
plot(fd,sec,xlim=c(-4.2,1.1),ylim=c(0.5,2.5),col="white",xaxt="n",xlab="",ylab="",yaxt="n")
segments(lb,sec,ub,sec,lwd=3,col="grey50")
points(fd,sec,pch=20,cex=2)
abline(v=0,lwd=2,lty=3)
axis(1,at=c(-4:3),cex.axis=1.4)
mtext("Marginal Effects: Binary Turnover Variable",3,outer=T,line=-2,cex=1.6)
axis(2,at=c(1,2),labels=c("Low Mil Spending","High Mil Spending"),cex.axis=1.4,las=2)
mtext("Percentage Point Change",1,outer=F,line=2.4,cex=1.4)


##########
## Figure A5
##########

  ## Quantities based on "Figure A5 Inputs" from "Schub_ISQ_LeaderTurnover_MainAnalysis.do"
  
 ## Left Panel: Democracies
turns<-c(0,1,2)
high<-c(0.043,0.053,0)
low<-c(0.028,0,0.008)

par(mar=c(3,4,2,2))
par(oma=c(2,1,0,0))
plot(turns,high,pch=19,col="grey50",cex=3.8,ylim=c(-0.002,0.06),yaxt="n",ylab="",xaxt="n",xlab="")
points(turns,low,pch=9,col="black",cex=3.8)
axis(2,at=seq(0,0.06,by=0.01),labels=c(0,1,2,3,4,5,6),cex.axis=1.4)
mtext("War Occurrence (%)",side=2,outer=F,line=2.6,cex=1.7)
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)
mtext("Democratic Targets",side=3,outer=F,line=0.4,cex=1.7)
text(1,0.049,"high military spending",col="grey35",cex=1.3)
text(1,0.006,"low military spending",col="black",cex=1.3)
text(0.12,0.04,"n=69",cex=1.1)
text(1.13,0.057,"n=57",cex=1.1)
text(1.84,0.0,"n=61",cex=1.1)

 ## Right Panel: Non-Democracies
turns<-c(0,1,2)
high<-c(0.02,0,0)
low<-c(0.005,0.008,0.004)

par(mar=c(3,4,2,2))
par(oma=c(2,1,0,0))
plot(turns,high,pch=19,col="grey50",cex=3.8,ylim=c(-0.002,0.06),yaxt="n",ylab="",xaxt="n",xlab="")
points(turns,low,pch=9,col="black",cex=3.8)
axis(2,at=seq(0,0.06,by=0.01),labels=c(0,1,2,3,4,5,6),cex.axis=1.4)
mtext("War Occurrence (%)",side=2,outer=F,line=2.6,cex=1.7)
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)
mtext("Non-Democratic Targets",side=3,outer=F,line=0.4,cex=1.7)
text(0.51,0.02,"high military spending",col="grey35",cex=1.3)
text(0.51,0.002,"low military spending",col="black",cex=1.3)
text(0.12,0.024,"n=304",cex=1.1)
text(1.13,0.003,"n=97",cex=1.1)
text(1.84,0.000,"n=66",cex=1.1)
 


##########
## Figure A6
##########

  ## Past Turnover Predicts Future Turnover
turn<-read.dta("PredictTurnover_Schub_ISQ_LeaderTurnover.dta")
   # solschdum=binary sols change in that year
   # sols10lagtri=three level coding for numbers of sols changes in preceding 10 years

  ## Left Panel: Descriptive
mat<-matrix(NA,nrow=3,ncol=3)
for(i in 1:3){
grab<-turn$solschdum[turn$sols10lagtri==(i-1)]	
mat[i,1]<-(i-1)	
mat[i,2]<-length(grab)
mat[i,3]<-mean(grab)	
}
means<-mat[,3] ## observational mean chisols
   #prior10=0 --> 4.7%; prior10=1 --> 11.4%; prior10>=1--> 14.1%
ct<-c(0:2)

par(mar=c(3,4,1,2))
par(oma=c(2,1,2,0))
plot(ct,means,pch=19,cex=2,xlim=c(-0.2,2.2),ylim=c(0,0.18),xlab="",xaxt="n",ylab="",yaxt="n")
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)
axis(2,at=seq(0,0.15,by=0.05),labels=c("0","5","10","15"),cex.axis=1.4)
mtext("Probability of Leader Turnover (%)",side=2,outer=F,line=2.6,cex=1.7)
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Descriptive Statistics",side=3,outer=T,cex=1.8,line=0)


  ## Right Panel: Modeled Predicted Probability
x<-glm(solschdum ~ sols10lagtri, turn,family="binomial")

annint<-cbind(turn,predict(x,turn,type="link",se=T))
annfull<-within(annint,{
	predprob<-plogis(fit)
	ll<-plogis(fit-(1.96*se.fit))
	ul<-plogis(fit+(1.96*se.fit))
	})

pred0<-mean(na.omit(annfull$predprob[annfull$sols10lagtri==0]))
pred1<-mean(annfull$predprob[annfull $sols10lagtri ==1])
pred2<-mean(annfull$predprob[annfull $sols10lagtri ==2])
ll0<-mean(annfull$ll[annfull $sols10lagtri ==0])
ll1<-mean(annfull$ll[annfull $sols10lagtri ==1])
ll2<-mean(annfull$ll[annfull $sols10lagtri ==2])
ul0<-mean(annfull$ul[annfull $sols10lagtri ==0])
ul1<-mean(annfull$ul[annfull $sols10lagtri ==1])
ul2<-mean(annfull$ul[annfull $sols10lagtri ==2])

main<-c(pred0,pred1,pred2)
lb<-c(ll0,ll1,ll2)
ub<-c(ul0,ul1,ul2)
ct<-c(0:2)

par(mar=c(3,4,1,2))
par(oma=c(2,1,2,0))
plot(ct,main,pch=19,cex=2,xlim=c(-0.2,2.2),ylim=c(0,0.18),xlab="",xaxt="n",ylab="",yaxt="n")
segments(ct,lb,ct,ub,lwd=3,col="grey40")
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Leader Turnovers in Prior Decade",side=1,outer=F,line=2.6,cex=1.7)
axis(2,at=seq(0,0.15,by=0.05),labels=c("0","5","10","15"),cex.axis=1.4)
mtext("Predicted Probability of Leader Turnover (%)",side=2,outer=F,line=2.6,cex=1.7)
axis(1,at=c(0,1,2),labels=c(0,1,"2 or more"),cex.axis=1.4)
mtext("Modeled Predicted Probability",side=3,outer=T,cex=1.8,line=0)



##########
## Figure A7
##########

  ## Peace costs decoupled from war costs
cost<-read.dta("PeaceWarCosts_Schub_ISQ_LeaderTurnover.dta")  
   ## Bilateral wars
    # milgdpestavg = average military spending before the war for the two belligerents
    # pcdeathBatsimple = average per capita battle deaths across the two belligerents 

plot(cost$milgdpestavg,cost$pcdeathBatsimple,pch=19,xlab="",ylab="")
abline(lm(pcdeathBatsimple ~ milgdpestavg,data=cost))
mtext("Mil Spending as % of GDP",side=1,line=2.5,outer=F,cex=1.5)
mtext("Battle Deaths per capita (000s)",side=2,line=2.5,outer=F,cex=1.5)
mtext("Bilateral Wars: Arms Spending and Fatalities",side=3,outer=T,line=0,cex=1.5)

  ## Correlation (-0.14)
t<-as.data.frame(cbind(cost$milgdpestavg,cost$pcdeathBatsimple))
to<-na.omit(t)
cor(to$V1,to$V2)











